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Abstract 

In this paper we propose several models that describe the dynamics of liquid films which are covered by a 
high concentration layer of insoluble surfactant. First, we briefly review the 'classical' hydrodynamic form 
of the coupled evolution equations for the film height and surfactant concentration that are well established 
for small concentrations. Then we re-formulate the basic model as a gradient dynamics based on an under- 
lying free energy functional that accounts for wettability and capillarity. Based on this re-formulation in the 
framework of nonequilibrium thermodynamics, we propose extensions of the basic hydrodynamic model 
that account for (i) nonlinear equations of state, (ii) surfactant-dependent wettability, (iii) surfactant phase 
transitions, and (iv) substrate-mediated condensation. In passing, we discuss important differences to most 
of the models found in the literature. 
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I. INTRODUCTION 



Small volumes of simple and complex fluids that occur naturally in biological contexts or that 
are employed in modern technology, such as e.g., in microfluidics, are often (partly) confined by 
a free surface that may be covered by surface active agents. These so-called surfactants may be 
tensids, lipids, certain nano-particles, or particular polymeric compounds. Because they decrease 
the surface tension of the free surface, gradients in their concentration correspond to gradients in 
the surface tension. These gradients result in tangential forces at the free surface that drive flows in 
the bulk liquid. This corresponds to the so-called solutal Marangoni effect, that is e.g., responsible 
for the tears of wine [73, 79]. 

All surface active agents are to some extent soluble in the bulk liquid, implying that a complete 
dynamical model needs to describe the motion of the bulk liquid, bulk concentration of surfac- 
tant, the surface concentration of surfactant and the adsorption/desorption processes that exchange 
surfactant molecules between the bulk liquid and the free surface. However, for many practically 
important surfactants, the bulk solubility is actually very small. Then one speaks of "insoluble 
surfactants" and only considers the dynamics of the surfactant that is adsorbed at the free surface. 
Here, we restrict our attention to such insoluble surfactants at concentrations at which no micelles 
are formed in the bulk liquid [6]. 

The governing transport equations that relate the material properties of the insoluble surfactant 
and the resulting hydrodynamic flow are well established for low values of the surfactant surface 
coverage F [19, 49, 59]. In this case, the linear equation of state 

7(r) =70 + 7rr (1) 

describes how the surface tension deviates from its reference value 7(0) = 70, for a bare free 
surface. The coefficient 7r is a material constant that is negative for most combinations of liquid 
and surfactant. The resulting tangential Marangoni force at the free surface is Vs7 = 7rVsr 
where Vs = (I — nn) ■ V is the derivative along the free surface and n is the unit normal vector. 
For any (linear or nonlinear) equation of state, the surface tension gradient Vs7 enters the tan- 
gential stress boundary condition of the momentum transport equation. The latter is accompanied 
by a transport equation for T that accounts for advective and diffusive transport of the surfactant 
[75, 84]. The resulting system of equations may be simplified in order to apply them to vari- 
ous physical situations such as, for example, the dynamics of surfactant-laden drops or bubbles 
immersed in (another) liquid [34, 47], free-standing soap films [24], liquid bridges covered by a 
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surfactant monolayer [51], surfactant-covered vertical falling liquid films [42], films on horizontal 
solid substrates [19, 40, 59], and drawn meniscii [72]. In particular, the latter geometry allows 
for an asymptotic treatment which results in a long-wave or lubrication description of the dynam- 
ics, via two coupled evolution equations for the film height and the surfactant surface coverage 
[19, 59]. In the following, we focus on this geometry, but we should emphasise that our main 
arguments also apply to the general case. 

Many works only treat the case of low surfactant surface coverages and employ the linear equa- 
tion of state in Eq. (1). The surface tension driven flow is then said to result from a linear solutal 
Marangoni effect. However, there is a growing literature where a similar approach is used to treat 
the dynamics of free surfaces covered by large concentrations of insoluble surfactants. It is com- 
mon practice to replace the linear equation of state (1) by a nonlinear one and leave all other terms 
in the dynamical equations unchanged. We argue below that this may result in governing equations 
that are thermodynamically inconsistent, since one must also amend the surfactant surface diffu- 
sion term. There should also be amendments to the basic equations which describe the influence 
of the surfactant coverage on wettability close to three phase contact lines of very thin films, the 
effects of phase transitions in the surfactant layer at high concentrations, and also the influence of 
a nearby solid substrate on such phase transitions that may result in substrate mediated condensa- 
tion (surfactant aggregation). The approach that we propose in this paper allows one to deal with 
all these cases in a thermodynamically consistent manner 

The structure of the paper is as follows: First, we review in section II the 'classical' thin-film 
hydrodynamic coupled equations of motion for a thin liquid film covered by a low concentration 
surfactant. Then, in a preparatory step we decouple and re-formulate the two individual equations 
in a "thermodynamic form". In particular, section III gives the gradient formulation of the evolu- 
tion equation for a thin-film of pure liquid on a solid substrate, for the case where capillarity and 
wettability are the dominant influences, while section IV briefly reviews the classical diffusion 
equation and places it in the thermodynamic context that we employ. In section V the full coupled 
system is re-formulated as a gradient dynamics based on an underlying free energy functional. 
This thermodynamic form is used in section VI to extend the thin-film model to consistently ac- 
count for (i) nonlinear equations of state, (ii) surfactant-dependent wettability, (iii) surfactant phase 
transitions, and (iv) substrate mediated condensation. We also note some differences to the models 
in the literature. Finally, section VIII concludes and discusses the hmitations of our approach. 
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FIG. 1 : Sketch of a surfactant covered liquid film. 



II. THIN-FILM EQUATION FOR LOW SURFACTANT SURFACE COVERAGE 

If a hydrodynamic system involves a free surface that is covered by an insoluble surfactant, 
the boundary conditions for the momentum equation have to be supplemented by an evolution 
equation for the surfactant concentration on the free surface that accounts for transport of the sur- 
factant by advection and diffusion and also for shape changes of the surface that act as effective 
source/sink terms [49, 75]. This equation must be solved in conjunction with the hydrodynamic 
equations and boundary conditions for the liquid film. These equations can be greatly simplified 
for the case of a thin film of liquid on a solid substrate. If all quantities in the film vary over 
distances with a length scale parallel to the substrate that is large as compared to all length scales 
perpendicular to it, one may make a long- wave approximation [19, 59] to obtain coupled evolution 
equations for the film thickness profile h(r, t) and the surfactant surface coverage profile r(r, t), 
which is a dimensionless surface packing fraction (or concentration), where r = (x, y) is a carte- 
sian coordinate over the surface. The surface concentration is defined as r(r, t) = Pp{r, t), where 
p(r,t) is the surface number density (number per area) and P is the surface area per surfactant 
molecule when the surfactant molecules are at maximum packing on the surface (i.e. / is a molec- 
ular length scale), so that close packing corresponds to F = 1. For the three-dimensional physical 
situation illustrated in Fig. 1 , the equation for the film height is 



dth 



-V- 



V (7oA/i -padd(/i)) 



- V- 



Irh' 
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Note that this equation is equally valid if one uses p instead of the dimensionless V or indeed any 
other measure of the surfactant surface density. Only the quantity 7r needs to be redefined so that 
the product of 7r with the surface density yields a quantity with the dimensions of an energy per 
area [c.f. Eq. (1)]. As a result, many papers in the literature do not mention what units they choose 
for r. The time evolution equation for V is 



ar = -V ■ 



V (7oA/l -padd(/i)) 



V- ( ^^^vr ) + V- (wr), (3) 



2r7 

where 70 is the liquid-gas surface tension and rj is the dynamic viscosity of the pure liquid. Partial 
derivatives with respect to time and space are denoted dt and dx, respectively, V = {dx,dy) is 
the planar gradient operator and A = dxx + dyy is the Laplace operator. The mobility Q{h) = 
h^/3r] results from Poiseuille flow in the film without slip at the substrate. The pressure p = 
— 7oA/;,+padd(^) contains the Laplace surface curvature contribution to the pressure and additional 
contributions such as a hydrostatic or a disjoining pressure [23, 24, 40, 55, 59]. Note that the latter 
is normally assumed to be independent of T. Exceptions are discussed below. The diffusive 
transport of the surfactant in Eq. (3) follows from Fick's law for the flux Jdis = —DVT. In most 
papers, it is assumed that the diffusion constant D does not depend on the surfactant concentration, 
i.e., the term V ■ (DVT) in Eq. (3) becomes DAT. 

Note that Eq. (3) is an equation obtained in the long-wave approximation and therefore does 
not include the source-like surface dilatation term. Different forms for such a term are discussed in 
Refs. [13, 62, 75]. For the same reason, the V operator in the diffusion term is the planar operator 
and not the operator Vs that acts tangentially to the free surface. To extend the present ideas to 
more general geometries, these contributions must be taken into account. 

To obtain Eqs. (2) and (3), we have related the surfactant surface coverage T to the surface 
tension 7 by the linear equation of state in Eq. (1), i.e., a linear solutal Marangoni effect is assumed. 
In deriving Eq. (2), one also assumes that 70 ^ 7r(ro — T) and that therefore the Laplace pressure 
term (— 70A/1) only depends on 70. 

To incorporate effects of high surfactant concentration, the equations are often extended by 
translating 7rVr back into V7(r), and then replacing the linear equation of state in Eq. (1) by 
some non-linear equation of state. There are problems with doing this, as we show below. An- 
other extension is to incorporate a surfactant-dependent wettability into the evolution equations 
(2) and (3). This is sometimes done in an ad-hoc manner by simply replacing Pa.dd{h) by some 
Padd(^, r). However, it turns out that this leaves the equation incomplete and may even result 



in qualitatively incorrect predictions. After re-formulating the evolution equations as a gradient 
dynamics in section V, we discuss such extensions in Section VI. First, however, we introduce 
the gradient formulation for the decoupled thin-fikn equation (section III) and surfactant surface 
diffusion equation (section IV). 



III. THIN-FILM OF PURE LIQUID - EVOLUTION EQUATION AS A GRADIENT DYNAMICS 



It was noted some time ago that the time evolution equation for the film thickness in the case 
without surfactant [Eq. (2) with F = 0] can be written in a variational form [57, 60]. This allows 
one to appreciate that Eq. (2) corresponds to a time evolution equation for a conserved order 
parameter field h{r, t) (cf. Ref. [48]) that follows a dissipative gradient dynamics governed by the 
following equation 



d,h 



V- 



SF' 

Jh 



(4) 



This equation describes how the field h evolves towards a minimum of the free energy functional 



F[h] 



dA 



(5) 



where f{h) = J Pa.dd{h)dh, the mobility function Q^h = h^/3r] (cf. also Ref. [76]) and dA is a 
cartesian area element along the substrate. 

Note that the free energy in Eq. (5) corresponds to the one that is obtained making a small 
slope approximation in the free energy F = J f{h)dA + J 'jodS, where the surface element 
dS = a/1 + (V/i)^ dA = ^dA is approximated using ,^ ~ 1 + (V/i)^/2. The constant part J ^odA 
of the free energy in Eq. (5) is normally omitted since it does not contribute to the dynamics, as 
one can see from Eq. (4). A similar formulation is given in the following section for the surfactant 
surface diffusion equation. 



IV. DIFFUSION EQUATION AS A GRADIENT DYNAMICS 

The diffusive transport of a species with small surface coverage F is described by the diffusion 
equation 

dtT = -V ■ Jdiff = DAT, (6) 

where the diffusive flux is given by Pick's law Jdifj = —DW. The time evolution equation for 
the surfactant density in Eq. (3), reduces to the diffusion equation in Eq. (6) in the limit when 
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the liquid film thickness h{r,t) is a constant and the coefficient 7r = 0; i.e. when there is no 
Marangoni effect. 

The form in Eq (6) can easily obscure the underlying thermodynamics, which can be seen when 
this equation is formulated as a gradient dynamics based on the Helmholtz free energy for an ideal 
gas (i.e. a system of non-interacting particles): 



where k is Boltzmann's constant and T is the temperature. The transport equation for T is of the 
same form as Eq. (4) and reads [4, 5, 27, 53, 54]: 



where the mobility Qrr = DT. Here, D is the molecular mobility related to the diffusion process 
and may in principle depend on all the independent variables, although in the following we will 
assume that D is constant. 

The equivalence of the formulation in Eq. (6) and in Eqs. (7) and (8) is easily established. 
The advantage of the gradient dynamics form in Eq. (8) is the "built-in" straightforward way 
to extend the description, e.g., to incorporate attractive forces between the diffusing molecules 
and the effects of higher concentrations. For instance, replacing the functional in Eq. (7) by the 
one discussed by Cahn and Hilliard [12] (their equation (2.4) with (3.1)) results in the non-linear 
diffusion equation (or Cahn-Hilliard equation), Eq. (9) of Ref. [10] (when Qrr = DT (1 — F) 
is expanded about F = 1/2 and only the lowest order term is kept). More recently, Marconi 
and Tarazona [53, 54] showed that one can derive Eq. (8), starting from over-damped stochastic 
equations of motion for the (surfactant) particles. They showed that the diffusive fluid dynamics is 
described by Eq. (8), taken together with a suitable approximation for the Helmholtz free energy 
functional taken from equilibrium density functional theory [27, 36]. This so called dynamical 
density functional theory [4, 5, 53, 54] is now a growing body of work, allowing one to go beyond 
Cahn-Hilliard theory, and to develop a theory which includes a microscopic (on the scale of the 
particles) description of the dynamics of particles suspended in a fluid medium, or in the present 
case of surfactant particles on the surface of the liquid. 

To our knowledge, no gradient dynamics formulation has yet been given for the evolution of 
a thin-film covered by insoluble surfactant as described by Eqs. (2) and (3). Since the system is 
relaxational, i.e., there is no energy influx, a variational formulation in terms of a pair of coupled 




(7) 



dtT = V- QrrV 



6F 
5f 



(8) 
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evolution equations for two conserved order parameter fields must exist and in fact is presented in 
the following section. 



V. EVOLUTION OF A SURFACTANT-COVERED FILM AS GRADIENT DYNAMICS 

To construct a gradient dynamics description for the full coupled system, Eqs. (2) and (3), we 
start by considering the Helmholtz free energy for the system F[h, T], that also turns out to be the 
Lyapunov functional for a surfactant covered thin liquid film. It contains contributions that result 
from wettability (adhesion), expressed in terms of the height profile of the film, J f(h)dA, and 
contributions from the surface J g{V)dS, where g{V) is the Helmholtz surface free energy density 
(i.e. an energy per area), and dS = a/1 + {'Vh)'^dA = ^dA ^ [1 + (Vhy /2]dA is a surface ele- 
ment. The contribution of the height profile are similar to those in Eq. (5). In the parameter regime 
where the linear equation of state (1) for the surfactant is valid, i.e. when the surfactant density 
is low, then the surfactant layer corresponds to a two-dimensional gas of surfactant molecules on 
the film surface. The corresponding contribution to the free energy density is the entropic (ideal- 
gas) term (kT/P) r(log T — 1) that on its own leads to a diffusion equation for F, as discussed in 
section IV. However, one may add other contributions that are relevant at higher concentrations, 
resulting from the interactions between surfactant molecules. Such contributions are discussed 
below in section VI. For the non-interacting case we obtain 



It should be noted that variations in h and in F are not independent: If locally the slope of h 
changes, the area of the liquid surface changes and thus the surface coverage F may change without 
any surfactant transport. To derive evolution equations that are of a gradient dynamics form, one 
needs a concentration variable that is independent of the film height profile h. We introduce a 
surface coverage F, as sketched in Fig. 2, that corresponds to the coverage F 'projected' onto the 
flat substrate (thus, in principle, F can become larger than one). It is given by 




(9) 



where 



^(F)= 7o + — F[log(F)-l]. 



(10) 



F 



dS^ 
dA 



F = ^F. 



(11) 
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FIG. 2: Sketch that indicates the relation between the projected surface coverage T and the coverage T on 
the modulated free surface, defined in Eq. (11). 

Using F[h, T/^] from Eq. (9), the long-wave hydrodynamic equations (2) and (3) are equivalent 
to the following general form for the time evolution equations 



8h ST 



dth = V ■ 
dtf = W- 

where the symmetric positive definite mobility matrix is 



^ ^5F ^ ^6F 
<3/^^V— + QrrV^ 
5h 5Y 



(12) 
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\Qhr 
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h^r 

2ri 



+ DT 



(13) 



Note that we have written Q in terms of T and h. The justification for using T and not F in the 
long-wave approximation mobilities will become clear below. To fully appreciate the equivalence 
of Eq. (12) [with Eqs. (9), (10) and (13)] to Eqs. (2) and (3), we calculate the variations of the free 
energy in Eq. (9): 

5F 



Sh 
5F_ 

5V 



dhf{h) - V ■ [uVh] 
9\ 



(14) 



where we have introduced the local surface grand potential density u = g — Tg'. Inserting these 



9 



results into Eq. (12), we obtain the following time-evolution equations: 

dth = V- 



'f^V[d,f{h)-V-iuVh)] + ^Vg' 



ar ^ ar = v ■ 



■V [dhfih) - V ■ {uVh)] + ( — + DV]Vg' 



(15) 



where we have used the fact that (V/i)^ ^ 1 to approximate ^ ^ 1 in the left-hand side of the 
second equation, as is appropriate in the long-wave limit. Used at this stage, this approximation 
leads to r ^ r. Note, however, that this approximation should be applied with caution [85]. The 
equations for dth and dtV exactly correspond to the hydrodynamic model [Eqs. (2) and (3)] if (i) 
one identifies the local surface grand potential density u with the surface tension 7 (see further 
discussion in Section V A below) and (ii) employs Eq. (10) for g(T), i.e., one only includes the 
ideal-gas contribution to the free energy, which is a reasonable approximation to make for low 
concentrations of surfactant. As a result 



V7 = -VV5F/5f = TVg' = {kT/f)VT, 



(16) 



and so one finds that the diffusion coefficient D and the solutal Marangoni coefficient 7r in Eqs. (2) 
and (3) are given hy D = kTD/P and 7r = —kT/P, respectively. 

Note that the hydrodynamic community often assumes that the change of the surface tension 
with concentration is small as compared to the reference surface tension 70 and therefore only uses 
7o in the Laplace pressure term, i.e., V ■ (uVh) ^ 7oA/i. 



A. Equation of state - surface tension 

Before presenting several extensions to the hydrodynamic equations (2) and (3), based on our 
reformulation in Eq. (15), we put our gradient dynamics formulation in its proper thermodynamic 
context. First, we discuss the equation of state that relates surface tension 7 to surfactant concen- 
tration r, and show that the thermodynamic and hydrodynamic approaches are fully consistent for 
any dependence g(T). 

To this end, we review first some elementary considerations concerning the relation between 
the surface tension and the surface free energy density. The latter is defined as the excess free en- 
ergy per unit area that is due to the presence of a surface. In analogy to bulk thermodynamics, this 
excess may be defined for different thermodynamic ensembles and may therefore depend on dif- 
ferent surface thermodynamic variables. Moreover, this surface excess acts as the thermodynamic 
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potential for surface variables (for an enlightening general discussion of surface excesses, see Ref. 
[11]). The surface tension is the derivative of this surface thermodynamic potential with respect 
to the area - it is the equivalent of the pressure in bulk thermodynamics (up to a sign convention: 
a positive pressure generates an outward force on the container walls, whereas a positive surface 
tension creates an inward force on the lines bordering a surface element). 

It is easy to show that the surface tension is always given by the surface excess grand potential 
density, regardless of whether the surface free energy is defined in the canonical or grand-canonical 
ensemble. In the canonical case, the surface free energy is the surface (excess) contribution to the 
Helmholtz free energy, Fsuri = Sg{T), for a surface element of area S. The variation of this free 
energy is dFsurf = g(r)dS + Sg'(T)dT. The second term arises from the fact that the variation in 
the canonical ensemble has to be taken for a fixed number of surfactant molecules, N = ST/P 
(recall that T = pP), and thus the variation of the surface area creates a variation in the local 
concentration equal to dT = d{PN/S) = -{I'^N/ S'^)dS = -(T/S)dS, and the surface tension 
becomes 7 = co{T) = dF^^^i/dS = g{T) — g'{T)T, where we have introduced the surface grand 
potential density uj{T) = g{T) — pT, with /i = g'{T) = dg/dV the chemical potential [86]. 
Alternatively, if the surface excess is defined in the grand-canonical ensemble, the surface free 
energy is directly given by Fgurf = Suj. Now, the variation has to be taken at constant chemical 
potential (the surface element is connected to a reservoir of surfactant molecules). Since F is a 
function of the chemical potential in the grand-canonical ensemble, it remains fixed, and therefore 
the surface tension is directly given by 7 = dFsuvi/dS = co [87]. 

For flat interfaces with a small surfactant concentration, the surface-related part of the local 
Helmholtz free energy g is given by Eq. (10), and so the chemical potential fi = (kT/P) logF. 
With this, the surface tension becomes 

j = oj = giT)-pT (17) 
kT 

= 70 - ^r, (18) 

i.e., one recovers the linear dependence, Eq. (1), used in hydrodynamics with 7r = —kT/t^ and 
Fq = 0. Note that if the surface tension is defined (incorrectly) as the local Helmholtz free energy 
g{V), the logarithmic terms entail that one does not recover the linear dependence in Eq. (1). 

By identifying the surface tension with the local grand potential, the thermodynamic and hydro- 
dynamic formulations are fully consistent for any convex local g{T). In the hydrodynamic formu- 
lation [Eqs. (2) and (3)] the Marangoni force contributes to the advective flux as —{h'^/27]) V7(F), 
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and normally the equation of state 7(r) is given directly. With •y = u = g — fiT, the Marangoni 
term becomes (h'^T/2r]) {drrd) VF. This expression is identical to the term one obtains in the 
variational formulation (12), i.e., Qrh^{SF/6f) = {h^T/2ri) V{6F/6f). For any F of the form 
in Eq. (9), this equals {h'^T/2r]) (drrd) VF. The equivalence of the Marangoni term in the hy- 
drodynamic and thermodynamic formulation is valid for any convex local giT). In other words, 
the surface tension gradient V7 may be expressed either as or as —Wg'. This implies that 
Eq. (15) may be written as 



dth 
dtT 



V 
V 



^V[dJ{h) - V ■ (7V/1)] - ^V7 

[dnf{h) - V ■ (7V/1)] - + D ) V7 



(19) 



Thus, the diffusion term is expressed in terms of V7 and the molecular mobility D. Note, however, 
that this argument no longer holds if the free energy functional contains non-local terms in F, such 
as, e.g., (VF)^. Then the formulation in Eq. (19) can not be used and one must start directly with 
Eqs. (12). 



VI. EXTENSIONS 



Up to this point, we have presented a gradient dynamics re-formulation of the hydrodynamic 
long-wave model for the evolution of a thin-film that is covered by a low concentration of insoluble 
surfactants, that has several connections to the approach taken in dynamical density functional 
theory, with a local approximation for the free energy. 

The present re-formulation really demonstrates its advantages when seeking to make (com- 
mon) extensions of the hydrodynamic model, such as to incorporate nonhnear equations of state, 
surfactant-dependent wettability, phase transitions at high surfactant concentrations, or substrate 
mediated surfactant condensation. Such effects are often included into the hydrodynamic formu- 
lation (2) and (3) in an ad hoc manner that may result in the omission of important terms and 
sometimes lead to qualitatively incorrect behaviour. The thermodynamic variational framework 
presented here, i.e., Eqs. (15), allows us to make extensions stemming from changes (extra terms) 
in the free energy functional (9) in a systematic and thermodynamically consistent manner. In the 
following, we discuss several examples. 
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A. Nonlinear equation of state 



The most common extension is to replace the linear equation of state in Eq. (1) by var- 
ious nonlinear expressions. Examples include the exponential (Gaussian) dependence 7 = 
Ci + C2exp(-rVC3) [74]; a smoothed step-wise change 7 = + C2tanh[C3(r - 1)] 
[35, 51, 52]; Langmuir-Szyszkowski (7 = Ci + C2log(l - T)) [34, 39, 47, 51, 61] or Frumkin 
(7 = Ci + C2log(l — r) + CsF^) [14, 51, 61] equations applied to insoluble surfactants; ex- 
pressions related to power laws, such as, e.g., 7 = Ci + C2(l + C^T)'^ [8, 32, 33, 40, 83]); and 
fits to experimentally obtained isotherms [9]. In all cases the C's represent various constants. In 
most works, the extension is done by solely replacing the parameter 7r in Eqs. (2) and (3) by the 
function 7' = d'j{T)/dT. This, however, does not take into account that in addition to the convec- 
tive transport due to the Marangoni force, the diffusive transport of the surfactant is also affected 
when the underlying free energy functional changes. Thus, most works assume that the diffusion 
constant D remains independent of the surfactant concentration, even when working with highly 
nonlinear equations of state. An exception is Ref. [8] that uses a D(T). 

Based on our thermodynamic reformulation, the proper relation between 7 and D (discussed at 
the end of section V A) results in the amended hydrodynamic equations (19) that are expressed in 
terms of 7(r). One may define a non-constant D{T) by enforcing Ficks's law 



to hold. From a comparison of Eqs. (3) and (15) one obtains DVT = DWg' = DTg"'VT, i.e., 
D = DVg" . With this D the hydrodynamic formulation (2) and (3) is consistent with the gradient 
dynamics form (15). To obtain D in terms of the equation of state, one differentiates the relation 
7 = — Vg' with respect to F. The result 7' = —Vg" implies 



Note, however, that in principle D itself might also be a function of F (and other state variables). 
The relations 7(F) and D(V) employed in Ref. [8] are only consistent with the general thermody- 
namic framework given here, if particular dependencies of D on F are assumed. Note that Fick's 
law in Eq. (20) is only true in the low density F — )► limit. More generally one should have 



(20) 



D = -Di. 



(21) 



Jdiff = -MV/i 



(22) 
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where M(T) is a mobility coefficient [4, 5, 22, 53, 54]. This form is universally valid, whereas 
Eq. (20) does not always hold (e.g. in the 'uphill diffusion' observed in spinodal decomposition). 

Note also that many authors correctly employ the hydrodynamic form in Eqs. (3) and (15) with 
a linear equation of state, but give the Marangoni flux in its general form as Jmar = (^^/S?]) V7(r) 
in combination with a surfactant-independent diffusion constant [19, 41, 56, 64, 82]. Whilst this 
approach is indeed correct for a linear equation of state, it should be stressed that it is not valid for 
arbitrary (nonlinear) equations of state. 



B. Surfactant-dependent wettability 

It is widely accepted that wettability depends on the surface density of the surfactants [7, 15]. 
However, the literature is less clear on how such effects may be incorporated in a hydrodynamic 
thin-film description by extending the model given in Eqs. (3) and (15). In all the contributions 
we are aware of, this is done by replacing the film-thickness dependent Derjaguin (or disjoining) 
pressure n(/i), that is contained in the pressure Padd(/i) in Eqs. (2) and (3), by a disjoining pressure 
n(/i, r) that depends on film thickness and surfactant concentration [18, 31, 37, 81]. The influence 
of surfactants on the various components of Derjaguin's pressure for thin-films are discussed in 
detail in Refs. [7, 21] in the context of free standing (soap) films. For a simple model for forces 
between surfaces with adsorbed layers see Ref. [38]. 

Based on our thermodynamic re-formulation, one can now see how the hydrodynamic equa- 
tions (2) and (3) must be amended to account for any dependency of the adhesion energy on 
film thickness and surfactant concentration. Replacing f(h) in the free energy functional (9) by 
f{h, r) results in the additional contributions drf to 5F/5T that affect both evolution equations. 
The resulting hydrodynamic form is 



dth = V ■ 
dtT = V ■ 



^v[d,f{h,r)-v-{^vh)] - ^[V7- rv9r/(/i,r)] 

V [d^J{K r) - V ■ (7V/1)] - + 5^ [V7 - TVdrfiK T)] 



(23) 



2r/ 

The extra terms should be interpreted as an additional contribution to the Marangoni force 
that must be taken into account for small film thicknesses. The effective Marangoni force is 
V7 + rV9r/(/i, r). It also becomes the effective driving force for diffusion of the surfactant. To 
our knowledge these terms have not been included in any of the thin-film evolution equations that 
model surfactant-covered ultrathin films. However, they are necessary in any model that involves 
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a surfactant-dependent Derjaguin pressure. Without them, the model may exhibit qualitatively 
incorrect behaviour, such as oscillatory instability modes [31], as discussed in more detail below. 
Furthermore, as the system evolves in time, it does not tend to the correct equilibrium state, par- 
ticularly in the contact line region. In short, to treat such effects correctly one must determine a 
suitable form for f(h, T), obtained from the surfactant-dependent Derjaguin pressures discussed 
in the literature [7, 15, 21]. 

C. Surfactant phase transitions 

A third example is the description of phase separating surfactant mixtures or phase transitions 
at high surfactant concentrations. The simplest case of surfactant molecules that shghtly attract 
each other is already addressed by the discussion above, as it only results in a nonlinear equation 
of state. For instance, for weakly attracting surfactant molecules, one replaces the purely entropic 
form in Eq. (10) by 5'(r) =70+7? r[log(r) — 1] — (a/2)r^, where the attraction strength parameter 
a > 0. This results in7 = 'jQ — (kT/P)T + (a/2)T'^ and so the effective diffusion constant depends 
linearly on F. 

The situation becomes more involved for surfactant layers that can undergo a phase transition 
when the concentration changes, e.g., between the gaseous and the liquid-expanded or between 
the liquid-expanded and the liquid-condensed phases [1, 71]. Beside a function g(T) that ac- 
counts for the particular surfactant isotherm, one also needs to incorporate a surface gradient term 
{k/2){V s^y^ in the free energy functional (9) to account for the finite width and line tension of 
the interface between the various surfactant phases. If a double- well potential is used for g, this 
amounts to a description of the surfactant layer using a convective Cahn-Hilliard-type equation. 
A similar approach is employed in Ref. [44] to describe a thin liquid film covered with an insol- 
uble surfactant in the vicinity of a first-order phase transition. However, as explained below, our 
formulation differs on a number of important points. 

As the algebra is involved, we illustrate this for the one-dimensional case, where x is the only 
spatial coordinate. The free energy functional is 



where we use ds = {1/C)dx, T = ^F, ds/dx = ^ and s is the arc-length coordinate along 
the free surface. Note that the final contribution to the integral in Eq. (24) is simply the term 
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(24) 



(/t/2) J{dsT)'^ds. We keep ^ = a/1 + {d^hy exact throughout the derivation and only use 
{d^hY ^ 1 at the end. Employing the approximation too early can lead to neglecting physi- 
cally essential terms, such as the Laplace pressure. Details of the calculation of the functional 
derivatives are given in Appendix A. The resulting expressions are 



SF 

Jh 
6F 

= g' - ndx^T (25) 

where we have used (dxhY ^ 1 [see appendix A - in particular Eqs. (A19) and (A20)] and 
CO = g — Tg'. The time evolution equations for h and T are obtained by substituting Eqs. (25) into 
Eqs. (12). On inspecting the resulting equations, one notices that we have again obtained the form 
in Eq. (19), but now the surface tension is 

7 = w = w-^(9,r)2 + fi:r9,,r. (26) 

Recall that above, for the case without the gradient terms in F, we had that d^^y = —Td^SF/ST 
[Eq. (16)]. It turns out that in the present case this result still holds. The surface grand potential 
density for the nonlocal case is uj = g + (K/2)((9^.r)^ — T5F/5T. This observation implies that 
with the proper definition of surface tension, the evolution equations in Eq. (19) are valid for both 
the extension to include nonlinear equations of state and the present extension that incorporates 
gradient terms in T in the free energy. 

This issue explains the differences between our formulation and that in Ref. [44], that starts 
from a hydrodynamic formulation somewhat similar to that in Eq. (19) [88]. We expect the for- 
mulation presented here to be useful for studying the dynamics of surfactant phase transitions on 
thin films, for the case of insoluble surfactants. For instance, incorporating gradient terms may 
enable one to explain the spatially non-monotonic distribution of a spreading surfactant drop that 
has been observed in recent experiments [29]. 



D. Substrate-mediated phase transitions 

As final example, we mention the so-called substrate-mediated phase transitions of surfactant 
layers that may occur when surfactant monolayers are transferred from a deep trough onto a solid 
substrate, i.e., during a Langmuir Blodgett transfer [50, 69]. Often, the substrate triggers a phase 
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transition from the liquid-expanded (LE) phase to the liquid-condensed (LC) phase of the sur- 
factant monolayer. Within the framework presented here, this transition may be described by 
replacing the surfactant contribution to the free energy (9), 5'(r), by a term that depends on both 
the surfactant concentration and film height, g(h,T). Doing this, one obtains additional contri- 
butions to the free energy variations: The g' has to be replaced by Org and a term dhg is added 
to 5F/5h in Eq. (14). For the full expressions, see appendix A. Such a surfactant concentration 
and film-height dependent contribution to the free energy is employed in Refs. [45, 46] to de- 
scribe substrate-mediated condensation, but without incorporating the additional dhg term. For 
their choice of g{h, F), the omission is of no major consequence; it only amounts to a redefinition 
of the parameters in the disjoining pressure. 



VII. CONSEQUENCES OF THE GRADIENT DYNAMICS FORMULATION 

The advantage of the gradient dynamics formulation, besides its thermodynamic consistency, 
is that one may readily use general results obtained for other systems having governing equations 
of the form of Eq. (12). Similar formulations exists, for instance, for two-layer thin-film systems 
(where the two conserved fields are the two film thicknesses) [67, 68] and for thin-films of so- 
lutions or suspensions (where the two conserved fields are the film thickness h and the effective 
solute layer thickness ^ = hcj), where </> is the vertically averaged concentration [77]). 



A. Lyapunov functional 



Just as in the above-mentioned cases, one can show that the free energy functional F[h,T] 
in Eq. (12) is a Lyapunov functional: The total time derivative of F[h,V] is dF/dt = 
/ (ff ~^ fp'^*"'^) Expressing the partial derivatives dth and dtV by the expressions in 
Eq. (12) and after integration by parts and assuming periodic or no-flux boundary conditions, one 
obtains 



dF 
~dt 



.6F 
Jh 



+ 2Q 



hT 



6h I 



6rJ 



Because [cf. Eq. (13)] 



det Q = + — — > 0, 



3r] 



V- 



.6£ 
6f 



dx. 



(27) 



(28) 
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and Qhh > and Qrr > 0, the quadratic form in Eq. (27) is positive definite and therefore 
dF/dt < 0, and F is a proper Lyapunov functional. Furthermore, one may identify the stationary 
solutions of Eqs. (12) with the extrema of F. 

B. Stability of flat films 

Next, we briefly discuss a general result for the linear stability of flat films h{x,t = 0) = 
ho, that are covered by a homogeneous layer of surfactant T{x,t = 0) = Tq. On a system of 
infinite size, one may decompose any fluctuation disturbances of the film height and surfactant 
concentration in this homogeneous state into Fourier modes and consider their time evolution. 
We employ the ansatz h{x, t) = /iq [1 + e exp {l3t + kx)] and T{x, t) = Fq [1 + ex exp + kx)] 
where k and (3(k) are the wave number and growth rate of the harmonic mode, respectively. The 
overall amplitude of the disturbance is e, while x is the amplitude ratio of the disturbances in the 
surfactant concentration and film thickness profiles. In short, the amplitudes may be written in 
vector notation as ex = e{ho, ^ro)^- 

Employing these ansatzes for h(x, t) and T{x, t) in Eqs. (12), and then linearising in e ^ 1, as 
is appropriate for small amplitude disturbances, leads to the following eigenvalue problem 

(J - /3I)x = 0, (29) 

where J is the non-symmetric Jacobian given by 

J = -A;2QoEo (30) 

and where Eg and Qq are the matrix of the second variations of F in Fourier space and the mobility 
matrix, respectively, both evaluated at Hq and Fq. Since det Q 7^ for /i, F > 0, Eq. (29) can be 
written as the generalised eigenvalue problem 

(A;2Eo + /3Qo')x = 0. (31) 

Because Eq and Qq ^ are both symmetric and Qq ^ is positive definite, one can deduce that all 
eigenvalues (3 are real [43], as one should expect for a variational problem. Inspecting Eq. (31) 
further indicates that the stability of the system is completely determined by the eigenvalues of 
Eq, i.e., by the second variations of the free energy functional. The stability threshold is given 
by detEo = 0. However, having Eq is not sufficient to obtain the actual growth rate of 
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the unstable modes (i.e. the dispersion relation) and the amplitude ratio x- These are obtained by 
solving Eq. (29). A remarkable effect that arises from the coupling of the two fields, i.e., when 
STSh 7^ 0, is that the system becomes unstable for a larger range of parameter values than the 
individual (decoupled) systems are on their own. This effect is discussed in many other contexts - 
see e.g. [16, 30, 67]. 

The fact that the present system has a non-diagonal mobility matrix Q in Eq. (13) (as do the 
models in Refs. [67, 77]) distinguishes it from many other systems with evolution equations for 
two coupled order parameter fields having a gradient dynamics. This means that the equations 
for both fields depend on variations of the free energy with respect to each of the fields h and 
r. Thus, for a non-diagonal Q and when 5'^F/5V5h ^ 0, the evolution of the two fields are 
coupled both through the free energy functional and through the dynamical mobility coefficients 
in Q. In contrast, many such systems have a diagonal mobility matrix, and then the evolution 
of the two fields is solely coupled by the off-diagonal term in the matrix of second derivatives, 
S^FjbVbh. Examples are the equations in Ref. [65], that describe the spinodal decomposition in 
ternary systems, the dewetting of nanoparticle suspensions in Ref. [70], the coupled demixing and 
dewetting of a binary mixture discussed in Ref. [17], the electric field driven surface instability 
of two air-gap separated polymer layers in a capacitor [2], and the model equations employed in 
Ref. [30] to describe the interplay between ordering and spinodal decomposition in binary systems. 

VIII. CONCLUSIONS 

We have proposed several amendments and extensions for models describing the dynamics of 
liquid films that are covered by high concentrations of insoluble surfactant. After briefly review- 
ing the 'classical' hydrodynamic form of the coupled evolution equations for the film height and 
the surfactant concentration profile, that are well established for small concentrations, we have 
re-formulated the model in three stages as a gradient dynamics. We refer to this as the "thermody- 
namic form" of the evolution equations. 

In the first stage, we have given the gradient dynamics form of the evolution equation for a 
thin film of a pure liquid on a flat substrate without surfactant, in the case where capillarity and 
wettability are the dominant influences. This formulation was discussed before, e.g. in Refs. [57, 
76], and is of the standard form suitable for conserved dynamics, of which a classic example is 
the Cahn-Hilliard equation [10] for the demixing dynamics of a binary mixture. In the second 
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stage, we have briefly reviewed the classical diffusion equation and have noted the existence of 
a gradient dynamics formulation that puts it in the context of nonequilibrium thermodynamics 
and dynamical density functional theory. Finally, in the third stage we have re-formulated the 
full coupled system of equations for the liquid film height and surfactant concentration profile 
in a gradient dynamics form, based on an underlying free energy functional that accounts for 
wettability, capillarity and entropic contributions for the surfactant. The resulting equations are 
equivalent to the hydrodynamic form for the case of a linear equation of state for the surfactant. 

Based on this thermodynamic re-formulation, we have proposed amendments to the basic hy- 
drodynamic model that account for four different physical effects that all may be included through 
changes to the free energy functional. In particular, we have extended the thin-film model to 
consistently account for (i) nonlinear equations of state, (ii) surfactant-dependent wettability, (iii) 
surfactant phase transitions, and (iv) substrate-mediated condensation. The ideas that we have 
presented can also be directly applied to films covered by monolayers of nano- sized particles that 
are not soluble in the liquid film [28], or any substance that remains on the surface of the liquid 
film. 

Our results indicate that nearly all long-wave models found in the literature that extend the hy- 
drodynamic equations for thin liquid films covered by insoluble surfactants by including non-linear 
equations of state are either not fully consistent or not complete. The most important differences 
between our model and those in the literature, as discussed above in section VI, are: 

(i) When incorporating a nonlinear equation of state, most authors fail to note that one must 
also amend the surfactant diffusion term in the governing dynamical equations. 

(ii) To account for a surfactant-dependent wettability, it is not sufficient to just adapt the Der- 
jaguin (or disjoining) pressure. The Marangoni and diffusion term must also be amended. 

(iii) To account for surfactant phase transitions, square gradient, or other non-local terms for the 
surfactant concentration must be incorporated into the free energy functional. In these, the 
gradient should be taken along the free surface. 

(iv) When incorporating terms to describe a surfactant phase transition that depend on the dis- 
tance between the film surface to the solid substrate (e.g., to describe substrate mediated 
condensation), the added coupling terms lead to additional terms in the equation of state 
and the Derjaguin pressure. 
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The corrections to models in the literature that result from following our approach will in many 
cases only result in quantitative (rather than qualitative) changes to the results, that may also be 
rather small. In some cases, however, the differences will be qualitative and significant. For in- 
stance, we believe that the oscillatory dewetting modes ("dewetting waves") described in Ref. [31] 
for one and two-layer films with surfactant, are present in the model as a consequence of a broken 
variational structure of the governing equations that stems from omitting terms in the equation 
for the time evolution of the surfactant concentration profile. Using the complete equations, all 
eigenvalues of the linearised problem are real, and thus all instability modes are monotonic. 

The various extensions that we have proposed may all be simultaneously included so as to 
account for more complex situations. The corresponding free energy functional is 



where f{h,T) is a generalised wetting interaction term and g(T,h) is a generalised local free 
energy of the surfactant on the free surface. Note that future work should identify the connections 
that must exist between these more general functions / and g because they both arise from the 
same molecular interactions; this is an issue that we have not touched upon here. 

Our approach may also be further extended to accommodate more general terms that one 
should expect to be present in the free energy functional, including non-local integral (convolu- 
tion) contributions to F, which are commonly used in dynamical density functional theory (DDFT) 
[4, 5, 53, 54, 66], which uses as input the free energy functionals coming from equilibrium density 
functional theory [27, 36]. 

We emphasise that points (i) and (iii) above are particularly important for a number of bio- 
physical systems such as, e.g., the description of the surfactant layers that reside on the aqueous 
thin-film of the lung lining [40], where (a) the equations of state that are used are strongly non- 
hnear and (b) experiments show that phase transitions frequently occur, e.g., in layers of porcine 
lung surfactant at the air-water interface at physiologically relevant concentrations and tempera- 
tures [58]. Similar results are found for calf lung surfactant, where an expanded-to-condensed 
phase transition is observed as the surfactant concentration is increased [25]. The dynamics of the 
ongoing surfactant phase transitions and their interaction with the hydrodynamics of the thin liq- 
uid film is highly important. The formulation and extensions we present here allows one to extend 
the 'classical' hydrodynamic thin-film models to include the more intricate thermodynamic effects 
based on equations of state (obtained from suitable free energy functionals) that are observed ex- 




(32) 
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perimentally. Note, however, that a real layer of lung lining is much more complicated than the 
idealised situations mentioned above [63], as it consists of mixtures of soluble surfactants. Work 
is currently underway to extend our approach to describe soluble surfactants [78]. The extension 
towards mixtures of surfactants is more straightforward, but including a third field in the system 
makes the algebra somewhat tedious. 

Further extensions that should be considered in the future concern the dynamical aspects. Here, 
we have employed the mobility matrix Q in Eq. (13), which is derived from the well known hy- 
drodynamic transport equations obtained for the simplest case of an insoluble surfactant with a 
hnear equation of state. Although we believe this approximation should hold over a large param- 
eter range, at very high surfactant concentrations one must make corrections. Although, slip at 
the solid substrate can easily be accounted for, it is not clear what changes to the mobility ma- 
trix Q should arise from incorporating surface viscosity effects and/or a no-slip condition at the 
surfactant-covered surface. 

Here we have only discussed the gradient dynamics formulation in the context of surfactant- 
covered liquid films on solid substrates, because the mathematical formulation is most convenient. 
However, it is important to note that most of the effects that we mention also occur in other 
geometries. A prominent example where our considerations also apply are soap films based on 
insoluble surfactants. This is important for many systems, such as those reviewed in [7], that 
involve, for instance, surfactant-dependent Derjaguin pressures and highly nonlinear equations of 
state. Another example are surfactant-covered drops of liquid immersed in another fluid where 
issue (i) is particularly relevant when, for example, studying the shear-driven deformation and/or 
breakup of such droplets [26, 34, 47, 61] or liquid bridges or threads [3, 20, 51, 80]. In this 
case, for instance, incorporating a nonlinear equation of state should also be accompanied by the 
corresponding amendment of the surfactant diffusion term. 

Appendix A: Variational calculus in the general case 

The free energy for the surfactant covered thin liquid film is 

F[hj/C\ = J |/(/,,f/o + ^7(f/^,/^)e + ^(v(r/o)'^| dA. (Ai) 

We define 

F[h, f = F^et + i^surf + i^grad (A2) 
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so as to be able to calculate the variations of the three terms in the free energy separately. In the 
following, we limit ourselves to the 1 -dimensional case and we often need to use the result: 

= yiTTW =-Ad.h) {d^^h) . (A3) 



1. Variations with respect to h 
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(A12) 
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2. Variations with respect to F 
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3. Collecting the terms 



The resulting expressions are 

1 



— = dhf + dhg - -r- 
on ax 
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e 



6F 



K 



^drf + drg - -^d^j: + -i{dj:){d^h)d^^h 



(A18) 



This seems the appropriate stage in the derivation to apply the long-wave approximation, i.e., to 

use {dxhY = e <^1. Therefore ^ 1 + {l/2)e'^ and one obtains 



-TT = dhf + dug - -r- 
on ax 



g - Tdrg - '- {d^vy + kTO^^T ) d,h - KT{d^hf{d^^h)d^T 



SF 
6f 



drf + drg - nd^x^ + K{dj:){dxh)dxxh 



(A19) 



The respective last term is 0{e'^) smaller than the other terms with prefactor k and can therefore 
safely be dropped, yielding 



— = dhf + dhg - -r- 
oh ax 



g - Tdrg - ^ {dxTf + kTOxxT ) d^h 



SF 

sf 



drf + drg - nd^xT 



(A20) 



Eqs. (25) in the main text are obtained by setting dhg = drf = 0, whereas the results in Eqs. (14) 
are obtained by setting dhg = drf = together with k = 0. 
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